Generic Conditions for Hydrodynamic Synchronization 
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Synchronization of actively oscillating organelles such as cilia and flagella facilitates self-propulsion 
of cells and pumping fluid in low Reynolds number environments. To understand the key mechanism 
behind synchronization induced by hydrodynamic interaction, we study a model of rigid-body rotors 
making fixed trajectories of arbitrary shape under driving forces that are arbitrary functions of 
the phases. For a wide class of geometries, we obtain the necessary and sufficient conditions for 
synchronization of a pair of rotors. We also find a novel synchronized pattern with a time-dependent 
phase shift. Our results shed light on the role of hydrodynamic interactions in biological systems, 
and could help in developing efficient mixing and transport strategies in microfluidic devices. 
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Introduction. The idea that hydrodynamic interac- 
tions at low Reynolds number can induce synchronization 
between active components with cyclic motion has been 
the subject of extensive studies since the pioneering work 
of G. I. Taylor [l|, and has culminated in a number of di- 
rect experimental demonstrations in recent years 
For example, the effective and recovery strokes of beat- 
ing cilia [5| are considered to be important for generat- 
ing their coordinated motion (metachrony) While 
resolving the intricate conformations of the elastic fil- 
aments is important for studying coordination in high 
density assemblies, it can be argued that at sufficiently 
low densities hydrodynamic interaction does not alter the 
beating pattern of the active filaments, such that they 
can be feasibly modeled as simple beads following fixed 
trajectories @, [l(| ■ The simplicity of this level of descrip- 
tion allows for complex many-body effects to be probed 
in large arrays of such beads with additional active inter- 
nal mechanisms [81 11, 121. 



There have been a number of recent studies on hydro- 
dynamic interaction between rotating or orbiting rigid 
bodies [j| [HI II|-15[. An emerging general picture sug- 
gests that rigid bodies making fixed trajectories do not 
easily synchronize. Rigid helices with parallel axes 1 141 ] or 
beads on circular trajectories [1(3] with constant driving 
torque do not synchronize, unless flexibility is introduced 
in the orientation of the rotation axis [TBI j or in the con- 
finement to the trajectory [l3[, respectively. Vilfan and 
Julicher [§] studied two beads on tilted elliptic trajecto- 
ries near a substrate, with a velocity-dependent driving 
force. They found that both the height-dependence of 
the drag coefficient and the eccentricity of the trajec- 
tories are necessary to stabilize the synchronized state. 
Ryskin and Lenz [lOj considered a more general model, in 
which each cilium is represented by a collection of beads 
connected to each other. Each bead makes a fixed tra- 
jectory of arbitrary shape under a driving force that is 
an arbitrary function of the phase. They applied the 



general framework to a variety of beating patterns that 
mimic the ciliary strokes, and found them to be able to 
stabilize traveling (metachronal) waves but not synchro- 
nized states. These results naturally lead to the follow- 
ing question: when do objects with fixed trajectories syn- 
chronize via hydrodynamic interaction? Here, we address 
this question by formulating generic and explicit criteria 
for hydrodynamic synchronization. 

We use a simple version of Ryskin-Lenz model in which 
each active object (rotor) is made of a single bead. We 
derive a necessary and sufficient condition for a pair of ro- 
tors to synchronize, in terms of the trajectory shape and 
force profile. We apply the obtained criterion to specific 
trajectories, and identify the form of the force profiles 
that cause synchronization. For circular trajectories, for 
example, we find the requirement that the logarithm of 
the force has a non-vanishing second-harmonic compo- 
nent of a specific sign, which originates from the second- 
rank tensorial nature of hydrodynamic interaction. We 
consider trajectories in the bulk and near a substrate, 
and those tilted relative to each other. We also develop 
an effective potential picture to examine the global sta- 
bility of the synchronized states, which reveals a novel 
synchronized pattern with a time-dependent phase shift. 

Dynamical equations. We consider a pair of rotors 
(indexed by i = 1,2) and assume that each is a spheri- 
cal bead of radius a that follows a fixed periodic trajec- 
tory Ti = Yi{4>i), where <pi = 4>i(t) is the phase variable 
with the period 27r [see Fig. [Ha)]. The bead is driven 
by an active force Fi = Fi(4>i) that is tangential to the 
orbit and is an arbitrary function of the phase. The hy- 
drodynamic drag force acting on the bead is given by 
gi = Cf v ( r i) — where ( — 6nr]a is the drag coeffi- 
cient [16]. and v(r) is the velocity field of the surrounding 
fluid. The tangential component of the drag force is bal- 
anced by the driving force acting on each rotor, namely, 
Fj +ti -gi = 0, where is the tangential unit vector of the 
orbit given by = r^/|r£| with r- = dxijd^i. Substitut- 
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FIG. 1: (a) A generic trajectory with its shape specified by 
R(</>). The bead is driven by the tangential force F(<f>). (b) 
Circular trajectories Qi ■ R(0) and Q2 ■ R(<£) with R(0) = 
6(cos (j>, sin (j>, 0) and the rotation matrices Qi , Q2 . Their ori- 
entations are specified by the unit vectors qi = Qi ■ e^. (c) 
Linear trajectories with R(</>) = R((f>)e x and their orienta- 
tions specified by qi = Qi ■ e x . The centers of the trajectories 
are both on the £-axis. 



ing the expression for the drag force with = into 
the force balance equation, we obtain the phase velocity 
as (pi = uji +U ■•v(r i )/\r' i \, where u)i((j>i) = i*i(^»*)/Cl r i I 
is the intrinsic phase velocity. The reaction force — gi 
exerted by the bead on the fluid generates the flow field 



v(r) = -^G(r,r 3 ).g^^CG(r,r 

j 



r'^, (I) 



where G(r, Tj) is the Oseen tensor describing the hydro- 
dynamic interaction in bulk fluid. On the RHS of Eq. 
(U}, we assumed |r — rj\ 3> a and retained the leading 
order term with respect to £G(r, r^) = 0(a/\r — Tj\). 
Substituting this into the above expression for the phase 
velocity, we arrive at the coupled phase oscillator equa- 
tion 
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(2) 



where = G(r i ,r :) ). 

We now assume that the two trajectories have the same 
shape but are oriented differently relative to the axis that 
connects their centers. We can write each trajectory as 
Ti{4>) = r.jo + Qi • R(</>), where r^o is the position of the 
center, R(</>) describes the shape of the trajectory, and 
Qi is a rotation matrix. We also assume that the cen- 
ter positions are on the cc-axis and are separated by the 
distance d (s> a) from each other, rio = (0,0,0) and 
r 20 = (d,0,0). Using r^) = |R'(0)|Q, • t(</>) with the 
unit vector t((j>) = R'(0)/|R'(0)| in Eq. we obtain 
the difference between the phase velocities as 



01-02 = W(0l) - w(</) 2 ) 



F{fh) f , s F(0i) f . , 



where 



w(0i) 



F(<h) 



#(<fc,<fc),(3) 



(4) 



is the intrinsic phase velocity, and we have introduced 
the coupling function 

H(fa,fa) = Qi-t(0 1 )-CGi 2 -Q 2 -t(0 2 ), (5) 

which is a dimensionless quantity of order 0(a/d). To 
examine the stability of the synchronized state, we set 
4>i = 4>(t) + 6(t), cf>2 = 4>{t) and linearize Eq. ((3]) with 
respect to the phase difference S(t), which gives the linear 
growth rate 



2^(0) 
F{4>) 



H{^4>). (6) 



Integrating the above result over the period T — 
J v d(f>/(f), we obtain the cycle-averaged growth rate as 
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T 



#[lnF(0)]'ff(0,0), 



(7) 



to the lowest order in the hydrodynamic coupling H . A 
stable synchronized state exists when T < 0. Equation 
(|7|) thus shows that a necessary condition for synchro- 
nization is that both the force profile F{<j>) and the hy- 
drodynamic coupling H (0, 4>) are not constant. For any 
given trajectory R(0) that gives a non-constant function 
H((j),(j>), we can prescribe a force profile F(4>) that sat- 
isfies the above condition. For example, the force profile 
F((f>) =F I + ^ dip (H(i>, i>)-TT) , with H being the 

cycle average of £?(</>, </>), makes T negative-definite to the 
leading order in the coupling. 

In order to calculate H (<fi, 0), we decompose the Oseen 
tensor into isotropic (I) and dyadic (D) parts as 



CGi 2 = Gi(r 12 )I + G D (r 12 



r 12 r 12 
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where G/(r) = Gc(r) = 3/4r and we have used 
r i2 = ri - r 2 = -de x + Qi • R(0i) - Q 2 • R(<fe) 



(8) 



(9) 



When the characteristic dimension b = max|R(0)| of 
the trajectory is much smaller than the distance, we 
can approximate the hydrodynamic interaction kernel as 
£Gi2 ~ Gr(d)l + GD{d)e x e x . Under this approximation, 
the coupling function ([5]) becomes 



H( 



G D (d) [ qi • t(0)] [q 2 • t{<t>)} + const, (10) 



where qi = Qi • e^. Note that the diagonal part of 
the hydrodynamic kernel gives a constant contribution 
to H{<f>, cj)) and hence drops off from the integral J7]). We 
now examine a number of cases in more details. 

Circular trajectories. As the first example, let us con- 
sider the circular trajectory [see Fig. QJb)] 



TL(4>) = &(cos 4>, sin (j>, 0) . 



(11) 
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FIG. 2: Examples of the force profiles that act to synchronize 
two beads on circular trajectories aligned on the x-axis. (a) 
F{4>) = F {1 ~ | sin(20)]. (b) F(</>) = F [l + \ sin (0 + f )] . 



For this trajectory, we have |R' (</>)[ = b and t(</>) = 
(— sin 0, cos cj), 0). First we consider mutually parallel tra- 
jectories with Qi = Q2 = I. In this case, we have = e x 



and H((f),(t>) = G D {d) sm z 



^Gn(d) cos(20) + const. 



Note that the factor cos 2cf> represents the second-rank 
tentorial nature of the hydrodynamic kernel. To stabi- 
lize the synchronized state, we can use the force profile 
F{4>) = F o [l-Asin(20)] (F > 0,0 < A < 1), which gives 
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A- 



-r—,^-0(A 2 ) < [see Fig. EJa) for illustration]. 
We can also use F(<j>) = F [l + B sin (</> + § )] (F > 

0, -1< B < I), which gives T = -^§^B 2 + 0(B 4 ) < 
[see Fig. GJb)]. In general, the synchronized state is lin- 
early stable if and only if the Fourier expansion of In F((j>) 
has a negative coefficient for sin 2(f). 

Next, we consider rotated circular trajectories. For 
each trajectory (i = 1,2), the rotation operator 
that acts on (JTTJ) is parameterized by the Euler an- 
gles (a,,ft,7i) as Q t = M z (-f i )M x (/3 i )M z (ai), where 
M x (0) and M 2 (9) are the matrices of rotation by an- 
gle 9 around the x and z-axis, respectively. It gives 
qi = (cos a, cos 7i — cos Qj cos ft sin 7* , — sin a, cos 7, — 
cos cos ft cos 7^, shift sin 7^). For a rotation in the xy- 
plane (ft = 7* = 0), we get if (0, 0) = — \Goid) cos(20 — 
Oil — 0^2) + const., and synchronization is induced by, for 
example, the force profile F((f>) = Fq[1 — Asm(2<j) — oi\ — 
"2)] (0 < A < 1). Next, a rotation in the j/z-plane (on — 
7i = 0) gives if (c/f>, 0) = — ^Gr>(d) cos ft cos ft cos(20) + 
const. Finally, circular trajectories that are vertical 
to the xy-plane (o^ = 0, ft = ?•), give if (0, 0) = 
— \Gjj(d) COS71 cos 72 cos(20) + const. All of these cases 
yield similar conditions for synchronization in terms of 
the second Fourier coefficients of hxF{(f>), as in the case 
of non-rotated circular trajectories discussed above. Note 
that the decay rate to a synchronized state is indepen- 
dent of the size of the trajectory b at the leading order 
for circular trajectories. 

Linear trajectories. Next we consider the linear tra- 
jectory [see Fig. [Tfc)] 



R(0) = R(ct>)e a 



(12) 



that gives t((/>) = sgn[R'(<p)]e x , which amounts to a 
constant contribution to the coupling function (|10p and 
hence neither stabilize nor destabilize the synchronized 
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FIG. 3: Examples of the effective potential V(A). (a) 
For the circular trajectory with F((j>) = Fo[l — |sin(2</>)], 
the unique stable solution is found at A = 0. (b) For the 
elliptic trajectory R(0) = 6(cos 4>, ^ sin cf>, 0) with F(<f>) = 
F [l - ^sin(20) + |sra(4<£)], 

we obtain the bistable solu- 
tions A = ±Ao with Ao — 0.297T, and a metastable solution 
at A = 7T. (c) The stable solution A = Ao in (b) gives a 
time-dependent phase shift S — S(t) in the original gauge. 



state. However, 0(b/d) corrections to the hydrodynamic 
kernel give non-constant contributions. Substituting Eq. 
(|9]) into (jHJ and ([5]), and retaining the first order term 
with respect to R{4>), we obtain the coupling function 



H^cjy) = -2R&) 
G D (d) 



Gj(d)(qi ■ q 2 )Px + G' D (d)qi x q 2x p x 



{qix<\2 • p + feqi • p) 



const, (13) 



where p = qi — q2- Note that the coupling is constant 
when qi = q2, because the distance between the two 
beads is constant when the two trajectories are paral- 
lel. When they are not parallel, the phase dependence 
is proportional to R(<fi). For example, for perpendicular 
trajectories with qi = e x , q2 = e y and the orbital profile 
R{(j>) = &cos0, the synchronized state is stabilized if and 
only if the Fourier expansion of In F(<f)) has a positive 
coefficient for sin <j>. 

Nonlinear stability analysis. We can analyze global 
stability of the synchronized state by a nonlinear evolu- 
tion equation for the phase difference. To derive it, first 
we reparameterize the trajectory by the new phase vari- 
able $ = $>(<f>) that makes the intrinsic phase velocity 
constant: $'(0) • F((f>)/(\R! ((/>){ = 2tt/T. In this gauge, 
the phase difference A = $1 - $ 2 = ^{4>i) - $(<fc) obeys 
[see Eq. ©] 



A = 



2tt 



jfog) _ F($ 2 + A) 
\) F($ 2 ) 



F($ 2 



if($ 2 + A,$ 2 ), (14) 



where F($ 2 ) = F{(j) 2 ) and if($i,$ 2 ) = H(<j>i,<j> 2 ). We 
now average Eq. (| 14[) over the period < t < T, as- 
suming that A on the RHS is constant over a cycle, 
which is justified to the leading order in the coupling 
if [Tjj. We thus obtain the evolution equation in the 
form of A = —V'(A) with an effective potential V(A). 
The potential is calculated and plotted in Fig. [3] for two 
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examples, which are both non-rotated (Q, = I) and in 
the far-field (d 3> b). For the circular trajectory, the 
coupling function H contains only the second harmonic 
as a function of cj>, which results in either the in-phase 
(A = 0) or anti-phase (A = tt) synchronization depend- 
ing on the force profile, as illustrated in Fig. G2[a). For a 
more complex trajectory like the ellipse, more than one 
stable and/or metastable solutions can be obtained, as 
shown in Fig. EJb). Note that the non-zero values of A 
generally means non-constant phase shift S = <f>i — (f>2 in 
the original gauge, as Fig. [3Jc) shows. 

Synchronization near a substrate. Finally, let us con- 
sider the case where the rotors are suspended at height 
h from a flat substrate (located at z = —h). In this case, 
the hydrodynamic coupling is expressed by the Blake ten- 
sor [18j, which takes into account the no-slip boundary 
condition on the substrate. For simplicity, we restrict 
ourselves to the case where the trajectories are confined 
in the xy-plane and are separated by a relatively large dis- 
tance, namely d ^> h,b. In this case, we can use Eq. ((HJ 
with Gi(r) =0,G D (r) = 9h 2 /r 3 . While all the results 
discussed above hold true within the above restriction, it 
is straightforward to extend the analysis to more general 
and complex geometries. 

Discussion. Our analysis shows that the requirement 
for hydrodynamic synchronization is nontrivial but not 
difficult to meet, and that a wide variety of beating pat- 
terns do induce synchronization. For the example of cir- 
cular trajectories, we found that the logarithm of the 
force should contain second-harmonic component of a 
specific sign. This can be met if the force profile has 
the second harmonic directly, or, via frequency doubling, 
if it has the first harmonic. However, force patterns that 
only have harmonics higher than two cannot synchro- 
nize. Dependences on the trajectory shape and geome- 
try can be summarized by representing the action of a 
rotor at far distance by force multipoles at a fixed po- 
sition. For circular orbits, T is independent of the size 
of the trajectory, which indicates that the dominant in- 
teraction comes from force monopoles when they are in 
bulk fluid, or dipoles near a substratefeach made of the 
force monopole and its mirror image |18l|). Linear oscilla- 
tors do not couple at the first order of the force-multipole 
expansion. Their leading order coupling comes from the 
monopole-dipole interaction for beads that are in bulk, 
or dipole-quadrupole interaction near a substrate. Linear 
oscillators are also special in the sense that they do not 
synchronize if they are parallel. This could be relevant 
to the synchronization of microswimmers [l9| , which can 
be minimally modeled by a linear configuration of three 
point forces [2(|. We have also derived a fully nonlinear 
evolution equation for the phase difference, which enables 
us to study the global stability of synchronized states 
with or without phase shifts. We note that the presence 



of stable phase shifts has been recently observed in ex- 
periments on the beating patterns of the flagella of C. 
reinhardtii [3]. It is not difficult to extend the present 
analysis to helices (as a model of flagella) or other ex- 
tended objects to make comparison with experiments. 

In conclusion, we have derived a generic and explicit 
criterion for the trajectory shape and force profile that 
stabilize synchronized states. The criterion could be 
helpful in understanding the collective behavior of ac- 
tive biological organelles, and designing active microflu- 
idic components that could be tuned in and out of syn- 
chronized states using mechanical signals communicated 
via hydrodynamic interactions. 
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